# Who Cares if You Vote? Partisan agreement and injunctive norms of voting
# Edward Fieldhouse (1), David Cutts (2), Jack Bailey (1)
# (1) The University of Manchester, (2) The University of Birmingham


# 1. Housekeeping ---------------------------------------------------------

# The following code requires that you have first opened the associated
# project file. If not, click the button on the top right and open it
# ("Open Project...").

# Note also that the Bayesian methods we use below rely on having installed
# the probabilistic programming language Stan. If you have yet to install
# Stan, please see https://mc-stan.org/users/interfaces/.

# Load packages

library(tidyverse)
library(tidybayes)
library(magrittr)
library(ggpubr)
library(jbmisc)   # devtools::install_github("jackobailey/jbmisc)
library(haven)
library(brms)
library(here)


# Load data from the British Election Study Internet Panel. Note: If you
# plan to use this data, please make sure that you download the most recent
# version available at https://www.britishelectionstudy.com.

bes <- read_sav(here("_data", "BES2017_W14_Panel_v0.4.sav.zip"))



# 2. Transform data (respondent-level) ------------------------------------

# Select variables

bes <- 
  bes %>% 
  select(
    # Respondent-level covariates
    id,
    wave2,
    age = Age,
    female = gender,
    edu = edlevelW1_W6,
    edu2 = educationW1_W6,
    married = marital,
    att = polAttentionW2,
    pid = partyIdW2,
    sq = partyIdSqueezeW2,
    str = partyIdStrengthW2,
    
    # Discussant-level covariates
    rel_1 = relationshipName1W2,
    rel_2 = relationshipName2W2,
    rel_3 = relationshipName3W2,
    app_1 = discussantApprovalVoteName1W2,
    app_2 = discussantApprovalVoteName2W2,
    app_3 = discussantApprovalVoteName3W2,
    vote_1 = discussantVoteName1W2,
    vote_2 = discussantVoteName2W2,
    vote_3 = discussantVoteName3W2,
    pol_1 = discussPolDaysD1W2,
    pol_2 = discussPolDaysD2W2,
    pol_3 = discussPolDaysD3W2,
    class_1 = discussantClassName1W2,
    class_2 = discussantClassName2W2,
    class_3 = discussantClassName3W2) %>%
  filter(wave2 == 1)


# Z-standardise age

bes$age_z <-
  bes$age %>% 
  rescale()


# Convert female to factor with male as reference

bes$female <-
  bes$female %>%
  as_factor()


# Add "other" and "unknown" categories to education variable

bes$edu <-
  bes$edu %>%
  as_factor() %>%
  factor(levels = c("No qualifications",
                    "Below GCSE",
                    "GCSE",
                    "A-level",
                    "Undergraduate",
                    "Postgrad",
                    "Other",
                    "Unknown"))


# Populate "other" and "unknown" education categories from educationW1_6 variable

bes$edu[bes$edu2 == 18] <- "Other"
bes$edu[is.na(bes$edu) == T] <- "Unknown"


# Marital status

bes$married <-
  bes$married %>%
  as_factor() %>%
  as_dummy(c("Married", "Living as married", "Civil Partnership")) %>%
  factor(labels = c("Other",
                    "Cohabiting"))


# Political attention: mark "Don't know" as NA and convert to numeric

bes$att <-
  bes$att %>%
  mark_na(9999) %>%
  as.numeric()


# PID: Combine PID and PID Squeeze

bes$pid[bes$pid %in% c(10, 9999)] <- bes$sq[bes$pid %in% c(10, 9999)]


# Party ID

bes$pid <-
  bes$pid %>%
  as_factor() %>%
  factor(labels = c("Conservative",
                    "Labour",
                    "Liberal Democrat",
                    "SNP",
                    "Plaid Cymru",
                    "UKIP",
                    "Green Party",
                    "BNP",
                    "Other",
                    "None",
                    "Don't know"))

bes$sq <- NULL


# PID Strength

bes$str <- 
  bes$str %>% 
  as_factor() %>% 
  mark_na("Don't know") %>% 
  fct_explicit_na("None/Don't know") %>% 
  fct_relevel("None/Don't know")


# Convert data from respondent-level to dyad-level

bes <- 
  bes %>% 
  data.frame() %>% 
  reshape(
    varying = names(bes)[11:25],
    sep = "_",
    ids = "id",
    direction = "long",
    timevar = "disc"
  )



# 3. Transform data (referent-level) --------------------------------------


# Dyadic relationship, reference = "Other"

bes$rel <- 
  bes$rel %>% 
  as_factor() %>%
  fct_relevel("Other")


# Dyadic approval

bes$app <- 
  bes$app %>% 
  as_dummy(1) %>% 
  factor(labels = c("Other", "Yes"))


# Dyadic partisanship

bes$vote <- 
  bes$vote %>% 
  as_factor() %>%
  factor(labels = c("Labour",
                    "Conservative",
                    "Liberal Democrat",
                    "SNP",
                    "Plaid Cymru",
                    "Green Party",
                    "UKIP",
                    "BNP",
                    "Other",
                    "None",
                    "Don't know"))

opt <- c("None", "Don't know")

bes <- 
  bes %>% 
  mutate(
    prt =
      case_when(
        pid == vote & !(pid %in% opt) & !(vote %in% opt) & pid != "Other" & vote != "Other" ~ "Shared Partisanship",
        pid != vote & !(pid %in% opt) & !(vote %in% opt) ~ "Opposing Partisanship",
        pid == "Other" & vote == "Other" ~ "Opposing Partisanship",
        pid %in% opt & !(vote %in% opt) ~ "Respondent Non-partisan",
        vote %in% opt & !(pid %in% opt) ~ "Discussant Non-partisan",
        pid %in% opt & vote %in% opt ~ "Shared Non-partisan"
      ) %>% 
      factor(levels = c("Shared Non-partisan",
                        "Shared Partisanship",
                        "Opposing Partisanship",
                        "Respondent Non-partisan",
                        "Discussant Non-partisan"))
  )


# Days discussing politics with discussant

bes$pol <- 
  bes$pol %>% 
  mark_na(9999) %>% 
  as.numeric()


# Discussant Middle class

bes$class <- 
  bes$class %>% 
  as_factor() %>% 
  as_dummy("Yes, middle class") %>% 
  factor(labels = c("Other", "Middle Class"))


# Remove NAs listwise

bes <-
  bes %>%
  na.omit()


# Compute number of discussants by ID

bes <- 
  bes %>% 
  add_count(id,
            name = "n_disc")



# 4. Fit model ------------------------------------------------------------

# Fit multilevel approval model with brms

m1 <- brm(formula = 
            app ~ 1 + rel + pol + class + prt + age_z + female + edu + married + att + str + n_disc + (1 | id),
          family = bernoulli(link = "logit"),
          prior = 
            prior(normal(0, 1.5), class = "Intercept") +
            prior(normal(0, 0.5), class = "b") +
            prior(exponential(2), class = "sd"),
          data = bes,
          iter = 2000,
          chains = 4,
          cores = 4,
          refresh = 5,
          control = list(adapt_delta = .99,
                         max_treedepth = 15),
          file = here("_output", "app1_w2"))


# Inspect model

plot(m1, ask = FALSE)
summary(m1)
bayes_R2(m1)



# 5. Create plots ---------------------------------------------------------

# Use brms to compute dyadic social and partisanship effects conditioning
# on mean and reference values

me <- conditional_effects(m1,
                          effects = c("rel", "prt"))


# Create relationship plot

rel_plot <-
  me$rel %>% 
  mutate(
    effect1__ =
      effect1__ %>% 
      factor(
        levels =
          c(
            "Other",
            "Your spouse or partner",
            "Another relative",
            "A friend",
            "A neighbour or co-worker"
          ),
        labels = 
          c(
            "Other",
            "Your spouse\nor partner",
            "Another\nrelative",
            "A friend",
            "A neighbour\nor co-worker"
          )
      )
  ) %>% 
  ggplot(
    aes(
      x = estimate__,
      xmin = lower__,
      xmax = upper__,
      y = effect1__
    )
  ) +
  geom_errorbarh(height = 0,
                 size = .5) +
  geom_point(size = 1) +
  geom_point(size = .3,
             colour = "white") +
  geom_text(data = me$rel %>% 
              mutate(
                effect1__ =
                  effect1__ %>% 
                  factor(
                    levels =
                      c(
                        "Other",
                        "Your spouse or partner",
                        "Another relative",
                        "A friend",
                        "A neighbour or co-worker"
                      ),
                    labels = 
                      c(
                        "Other",
                        "Your spouse\nor partner",
                        "Another\nrelative",
                        "A friend",
                        "A neighbour\nor co-worker"
                      )
                  )
              ),
            aes(label = scales::percent(estimate__,
                                        accuracy = 0.1,
                                        suffix = ""),
                x = estimate__ + .002,
                y = effect1__),
            family = "Cabin",
            size = 2,
            vjust = 2,
            hjust = 0.7,
            inherit.aes = FALSE) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, .35)) +
  labs(title = "Dyadic Relationship",
       y = "",
       x = "Prob. of Discussant Approval") +
  theme_bailey() +
  theme(legend.position = "bottom",
        legend.margin = margin(t = 0, b = 0, unit = "cm"),
        legend.title = element_blank(),
        legend.text = element_text(size = rel(1)),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(hjust= .5))
rel_plot


# Create partisanship plot

prt_plot <-
  me$prt %>% 
  mutate(
    effect1__ =
      effect1__ %>% 
      factor(
        levels =
          c(
            "Shared Non-partisan",
            "Shared Partisanship",
            "Opposing Partisanship",
            "Respondent Non-partisan",
            "Discussant Non-partisan"
          ),
        labels = 
          c(
            "Shared\nNon-partisan",
            "Shared\nPartisanship",
            "Opposing\nPartisanship",
            "Respondent\nNon-partisan",
            "Discussant\nNon-partisan"
          )
      )
  ) %>% 
  ggplot(
    aes(
      x = estimate__,
      xmin = lower__,
      xmax = upper__,
      y = effect1__
    )
  ) +
  geom_errorbarh(height = 0,
                 size = .5) +
  geom_point(size = 1) +
  geom_point(size = .3,
             colour = "white") +
  geom_text(data = me$prt %>% 
              mutate(
                effect1__ =
                  effect1__ %>% 
                  factor(
                    levels =
                      c(
                        "Shared Non-partisan",
                        "Shared Partisanship",
                        "Opposing Partisanship",
                        "Respondent Non-partisan",
                        "Discussant Non-partisan"
                      ),
                    labels = 
                      c(
                        "Shared\nNon-partisan",
                        "Shared\nPartisanship",
                        "Opposing\nPartisanship",
                        "Respondent\nNon-partisan",
                        "Discussant\nNon-partisan"
                      )
                  )
              ),
            aes(label = scales::percent(estimate__,
                                        accuracy = 0.1,
                                        suffix = ""),
                x = estimate__ + .002,
                y = effect1__),
            family = "Cabin",
            size = 2,
            vjust = 2,
            hjust = 0.7,
            inherit.aes = FALSE) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, .35)) +
  labs(title = "Dyadic Party ID",
       y = "",
       x = "Prob. of Discussant Approval") +
  theme_bailey() +
  theme(legend.position = "bottom",
        legend.margin = margin(t = 0, b = 0, unit = "cm"),
        legend.title = element_blank(),
        legend.text = element_text(size = rel(1)),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(hjust= .5))
prt_plot


# Create plot of difference distribution

diff <-
  fitted(
    m1,
    newdata =
      tibble(
        id = NA,
        prt = c("Shared Partisanship", "Opposing Partisanship"),
        app = "Other",
        rel = "Other",
        pol = mean(bes$pol),
        class = "Other",
        age_z = mean(bes$age_z),
        female = "Male",
        edu = "No qualifications",
        married = "Other",
        att = mean(bes$att),
        str = "None/Don't know",
        n_disc = mean(bes$n_disc)
      ),
    summary = F,
    re_formula = NA
  ) %>% 
  data.frame() %>% 
  rename(
    shared = 1,
    opposing = 2
  ) %>% 
  mutate(diff = shared - opposing) %>% 
  pivot_longer(
    cols = everything()
  ) %>% 
  mutate(
    name = 
      factor(name,
                  levels = rev(c("shared", "opposing", "diff")),
             labels = rev(c("Shared\nPartisanship", "Opposing\nPartisanship", "Difference\nDistribution")))
  ) %>% 
  group_by(name) %>% 
  summarise(
    est = median(value),
    lower = quantile(value, 0.025),
    upper = quantile(value, 0.975)
  )

diff_plot <- 
  diff %>% 
  ggplot(
    aes(
      y = name,
      x = est,
      xmin = lower,
      xmax = upper
    )
  ) +
  geom_errorbarh(height = 0,
                 size = .5) +
  geom_point(size = 1) +
  geom_point(size = .3,
             colour = "white") +
  geom_text(data = diff,
            aes(label = scales::percent(est,
                                        accuracy = 0.1,
                                        suffix = ""),
                x = est + .002,
                y = name),
            family = "Cabin",
            size = 2,
            vjust = 2,
            hjust = 0.7,
            inherit.aes = FALSE) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1),
                     limits = c(0, .35)) +
  labs(title = "Difference",
       y = "",
       x = "Prob. of Discussant Approval") +
  theme_bailey() +
  theme(legend.position = "bottom",
        legend.margin = margin(t = 0, b = 0, unit = "cm"),
        legend.title = element_blank(),
        legend.text = element_text(size = rel(1)),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(hjust= .5))
  

# Group plots together

plot1 <- ggarrange(rel_plot, prt_plot, diff_plot, nrow = 1)


# Now we can render the plot as a jpeg file

jpeg(file = here("_output", "w2_approval.jpeg"),
     res = 300,
     width = 5.97,
     height = 2.2,
     units = "in")
plot1
dev.off()
